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Abstract 

In contrast to the original Kohn-Sham (KS) formalism, we propose a density functional theory 
(DFT) with fractional orbital occupations for the study of ground states of many-electron systems, 
wherein strong static correlation is shown to be described. Even at the simplest level represented 
by the local density approximation (LDA), our resulting DFT-LDA is shown to improve upon KS- 
LDA for multi-reference systems, such as dissociation of H2 and N2, and twisted ethylene, while 
performing similarly to KS-LDA for single-reference systems, such as reaction energies and equilib- 
£>~V rium geometries. Because of its computational efficiency (similar to KS-LDA), this DFT-LDA is 

Q-r applied to the study of the singlet-triplet energy gaps (ST gaps) of acenes, which are "challenging 

problems" for conventional electronic structure methods due to the presence of strong static corre- 
lation effects. Our calculated ST gaps are in good agreement with the existing experimental and 
high-level ab initio data. The ST gaps are shown to decrease monotonically with the increase of 
chain length, and become vanishingly small (within 0.1 kcal/mol) in the limit of an infinitely large 
polyacene. In addition, based on our calculated active orbital occupation numbers, the ground 
states for large acenes are shown to be polyradical singlets. 
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I. INTRODUCTION 



As the problem of solving the iV-electron Schrodinger equation quickly becomes in- 
tractable as the size of a system increases, the development of efficient and accurate electronic 
structure methods for large systems, continues being the subject of intense current inter- 
est. Over the past two decades, Kohn-Sham density functional theory (KS-DFT) [1,2] has 
become one of the most popular theoretical approaches for calculations of electronic prop- 
erties of large ground-state systems (up to a few thousand electrons), due to its favorable 
balance between cost and performance ^3 T | - Recently, its time-dependent extension, time- 
dependent density functional theory (TDDFT), has also been actively developed for treating 
electron dynamics and excited states of large systems with considerable success js, 9|. 

Although the exact exchange-correlation (XC) functional E xc [p] in KS-DFT has not been 
known, functionals based on the standard density functional approximations (DFAs), such 
as the local density approximation (LDA) and generalized gradient approximations (GGAs), 
can accurately describe short-range XC effects (due to the accurate treatment of on-top hole 
density), and are computationally favorable for large systems j^-7|. Although KS-DFAs have 
been successful in many applications, due to the lack of accurate treatment of nonlocality 
of XC hole, KS-DFAs can exhibit the following three types of qualitative failures: (i) self- 
interaction error (SIE), (ii) noncovalent interaction error (NCIE), and (iii) static correlation 
error (SCE). In situations where these failures occur, KS-DFAs can produce erroneous results 

n 

[10]. Therefore, resolving the qualitative failures of KS-DFAs at a reasonable computational 
cost seems to be the first step toward finding an efficient and accurate electronic structure 
method for large systems. 

The SIEs of KS-DFAs lead to drastic failures for problems such as barrier heights of 
chemical reactions, band gaps of solids, and dissociation of symmetric radical cations 
In TDDFT, SIE causes failures for problems such as Rydberg excitations in molecules and 
long-range charge-transfer excitations between two well-separated molecules [12| . The SIEs 
of KS-DFAs can be greatly reduced by hybrid DFT methods [13], incorporating some of the 
exact Hartree-Fock (HF) exchange into the KS-DFAs. Over the past twenty years, several 
hybrid functionals, such as global hybrid functionals fl^l . and long-range corrected (LC) 



iiO- 



hybrid functionals 
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20|, have been developed to improve the accuracy of E xc [p], extending 



the applicability of KS-DFT to a wide range of systems. 



The proper treatment of noncovalent interactions requires the accurate description of 
dynamical correlation effects at medium and long ranges, which cannot be properly cap- 
tured by KS-DFAs |2l| . In particular, for dispersion-dominated (van der Waals (vdW)) 
interactions, KS-LDA tends to overestimate the binding energies, while KS-GGAs tend to 
give insufficient binding or even unbound results. The NCIEs of KS-DFAs can be efficie ntly 



reduced by the DFT-D (KS-DFT with empirical dispersion corrections) schemes 22j-|25j, 



which have shown generally satisfactory performance on a large set of noncovalent systems 
The dispersion corrections can also be computed less empirically by the exchange- 
lole dipole moment (XDM) method 28] or by the local response dispersion (LRD) method 
29)]. Alternatively, a fully nonlocal density functional for vdW interactions (vdW-DF) 
can also be adopted to reduce the NCIEs of KS-DFAs. Currently, perhaps the most suc- 
cessful approach to taking into account nonlocal dynamical correlation effects is provided 
by the double-hybrid (DH) methods 3]]- 33], which mix some of the HF exchange and some 
of the nonlocal orbital correlation energy from the second-order M0ller-Plesset perturbation 
(MP2) theory {34 1 into the KS-DFAs. DH functionals have shown an overall satisfactory 
accuracy for thermochemistry, kinetics, and noncovalent interactions. In addition, the sharp 
increase in HF exchange in typical DH functionals has also greatly reduced the SIEs relative 
to KS-DFAs and conventional hybrid functionals. 

Systems with strong static (nondynamical) correlation effects, such as bond-breaking 
reactions, diradicals, conjugated polymers, magnetic materials, and transition-metal com- 
pounds, belong to the class of strongly correlated (SC) systems (multi-reference systems). 
Such a system is usually characterized by a small (or vanishing) energy gap between the 
highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital 
(LUMO), the HOMO-LUMO gap. Despite their computational efficiency, the accurate treat- 
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ment of SC systems poses a great challenge to KS-DFAs and hybrid DFT methods 
The DH methods may also be inadequate for SC systems, as the second-order perturbation 
energy components diverge to minus infinity for systems with vanishing HOMO-LUMO gaps. 
Within the framework of KS-DFT, fully nonlocal XC functionals, such as those based on 
random phase approximation (RPA), are believed to be essential for the accurate treatment 
of SC systems 6|, [7|. However, compared with KS-DFAs, hybrid functionals, and DH func- 
tionals, RPA-type functionals are computationally very demanding for large systems, and 
their applications to SC systems are too scarce to make a firm judgment on their accuracy. 



Recently, the SIEs of RPA-type 
one-electron system such as 



unctionals have been shown to be severe, even for a simple 



36|- 



Aiming to reduce the SCEs of KS-DFAs without extra computational cost, we propose 
a DFT with fractional orbital occupations, rather than a fully nonlocal XC functional in 
KS-DFT. The rest of this paper is organized as follows. In section II, we briefly describe 
the rationale for fractional orbital occupations. In section III, we describe the formulation 
of this DFT and explain how strong static correlation is described by this DFT. At the 
simple LDA level, the performance of our resulting DFT-LDA is compared with KS-LDA 
for several single- and multi-reference systems in section IV. Based on physical arguments 
and numerical investigations, the optimal parameter for this DFT-LDA is defined in section 
V. Our conclusions are given in section VI. 



II. RATIONALE FOR FRACTIONAL ORBITAL OCCUPATIONS 

For the exact DFT, the exact ground-state energy can be obtained, only when the exact 
ground-state density is available to insert into the exact ground-state energy functional. 
Therefore, the development of a generally accurate DFT method should involve not only 
an accurate approximation for the ground-state energy functional, but also an appropriate 
representation (possibly in terms of orbitals and occupation numbers) of the ground-state 
density. However, much less attention has been paid to the latter (representation of ground- 
state density) than to the former (ground-state energy functional). Due to the search over 
a restricted domain of densities, the exact ground-state density of interest may not be 
obtained within the framework of KS-DFT, in which case even the exact KS-DFT will fail 
i.e. the exact XC functional may not be different iable at the exact ground-state density) 
7fl. Noticeably, some of these situations occur for systems with vanishingly small HOMO- 
LUMO gaps (SC systems), indicating the close relationship between strong static correlation 
effects and representations of the ground-state density. Therefore, to accurately describe SC 
systems, it seems intuitive to focus on devising an appropriate representation for the exact 
ground-state density and a DFT associated with such a representation. 

A ground-state density p(r) is said to be interacting f-representable, if it can be obtained 
from a ground-state wave function of an interacting A-electron Hamiltonian for some ex- 
ternal potential v (r). The exact p(r) can be obtained by the full configuration interaction 



(FCI) method at the complete basis set limit (i.e. interacting f-representable) 37j, and can 
be compactly expressed in terms of the natural orbitals (NOs) {Xi( r )} an d natural orbital 
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occupation numbers (NOONs) {rii} 

oo 

PFCi(r) = ^ni|xi(r)| 2 , 
i=i 

where {rii}, obeying the following two conditions, 

oo 

^2m = N, < rii < 1, 



(1) 



(2) 



i=l 



are related to the variationally determined coefficients of the FCI expansion. As shown in 
Eq. ([T]), the exact p(r) can be represented by orbitals and occupation numbers, highlighting 
the importance of ensemble-represent able densities. 

By contrast, in KS-DFT |2j, p(r) is assumed to be noninteracting pure-state (NI-PS) 
f s -representable, as it belongs to a one-determinantal ground-state wave function of a non- 
interacting iV-electron Hamiltonian (KS Hamiltonian) for some local potential v s (r) (KS 
potential) 39M4ll|. Correspondingly, the KS orbital occupation numbers should be either 
or 1. As the Aufbau principle (i.e. filling the KS orbitals in order of increasing energy) 
should be obeyed, p(r) can be expressed as the sum of the densities of the lowest N occupied 
KS orbitals {&(r)} |a|: 

N 

Ei^wi 2 - ( 3 ) 



Pks-dftvJ 



i=l 



As ground-state densities of most nondegenerate atomic and molecular systems (e.g. closed- 
shell systems with sizable HOMO-LUMO gaps) are likely to be NI-PS unrepresentable, 
the commonly used XC functionals in KS-DFT are reliably accurate for these systems (as- 
suming that the SIEs and NCIEs of these functionals are not severe). However, if a one- 
determinantal ground-state wave function is insufficient to represent p(r), these XC func- 
tionals may not be reliably accurate, as they are all developed based on general properties 
of systems where the KS wave function is a one-determinantal wave function 



As shown by several researchers 
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44J ] . there are some reasonable ground-state densi- 



ties that are not NI-PS t> s -representable. Clearly, such densities cannot be obtained within 
the framework of KS-DFT. To remedy this situation, KS-DFT has been extended to en- 
semble DFT (E-DFT), wherein p(r) is assumed to be noninteracting ensemble (NI-E) -un- 
representable, as it belongs to an ensemble of pure determinantal states of the noninteracting 



KS system 45|, |46|. Correspondingly, p(r) can be expressed as 

oo 

Pe-dftW = ^^|0i(r)| 2 , (4) 
i=i 

where the occupation numbers {^}, obeying the following two conditions, 

oo 

5> = iV, 0<9i<l, (5) 

i=l 

are given by 

1, for 6i < €f 

9i= < Xj, for £» = e F (6) 
0, for €i > ep. 

Here is the orbital energy of the i th KS orbital <£ f (r), tp is the Fermi energy, and is 
a fractional number (between and 1). As can be seen, fractional orbital occupations can 
occur only for the orbitals at the Fermi level. 

In 1998, the close relationship between strong static correlation effects and non-NI-PS 
f s -representable densities was observed by Baerends and co-workers 42j, who studied the 
X E+ ground state of the C2 molecule (a system where the ground-state wave function is 
nondegenerate but has a strong multi- reference character), and investigated the possibility 
that the ground-state density p(r) may be NI-E v ^-representable (as assumed in E-DFT), 
rather than NI-PS t> s -representable (as assumed in KS-DFT). In their study, p(r) was first 
obtained from the highly accurate ab initio CI wave function method, wherein p(r) can be 
represented by the NOs {xi{r)} and NOONs {rii} in Eq. ([T|). An iterative method for the 
construction of the KS orbitals and the KS potential from the CI density was then adopted 
|47l |. combined with the constraint of integer occupations of the KS orbitals. The p(r) of 
C2 was found to be NI-PS t> s -representable at a bond distance shorter than the equilibrium 
distance, while being non-NI-PS t> s -representable at the longer bond distances, leading to 
non-Aufbau solutions with unoccupied KS orbitals having energies lower than those of the 
highest occupied KS orbitals (i.e. holes below the Fermi level). On the other hand, when the 
p(r) of C2 was represented by the ensemble solution during the above iterative procedure, 
no holes below the Fermi level were found, and the corresponding KS orbitals were close to 
the NOs. Based on their results, Baerends and co-workers argued that the KS orbitals 
(generated from Aufbau solutions) for a NI-PS t> s -representable ground-state density are 



comparable to the NOs , while the KS orbitals (generated from non-Aufbau solutions) for a 
non-NI-PS i> s -representable ground-state density can be distinctly different from the NOs in 
order to incorporate the effect of the configuration mixing on the ground-state density. They 
concluded that the ground-state density of a system with strong static correlation effects 
is likely to be non-NI-PS t> s -representable, in which case an ensemble representation (via 
fractional orbital occupations) of the ground-state density is crucial. Arguments in support 



of this are also available elsewhere 
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The idea of simulating strong static correlation effects by fractional occupation numbers 



(FONs) in DFT has spurred the development of the D 
restricted ensemble-referenced KS (REKS) method 53 



(FS-DFT) method 



3. 



T-FON method 



48 



52|, the spin- 



and the fractional-spin DFT 



35j . with great success for some SC systems. The practical imple- 
mentation of E-DFT and related methods has, however, been impeded by a few factors as 
follows. First, a double-counting of correlation effects may occur, when the conventional 
approximate XC functional is evaluated with an ensemble density in Eq. (Pflh rather than 
a NI-PS f s -representable density in Eq. ([3]). By taking into account the double-counting 
effects, the XC functional in E-DFT may need to be re-derived. Second, the computational 
cost of E-DFT is more expensive than that of KS-DFT, which makes E-DFT less practical for 
the study of large ground-state systems. Third, the computational cost of analytical nuclear 
gradients (if available) for E-DFT is more expensive than that for KS-DFT, which makes 
it a formidable computational task to perform geometry optimization of large molecules for 
E-DFT. 

In view of the above difficulties, we focus on the representation of the ground-state 
density from the exact theory in Eq. (CQ). Although the exact orbital occupation numbers for 
interacting electrons are intractable for large systems due to the exponential complexity, the 
approximate ones can, however, be properly simulated. Based on the statistical properties of 
strongly correlated eigenstates, Flambaum et al. argued that the distribution of occupation 
numbers (the microcanonical averaging of the occupation numbers) for a finite number of 
interacting Fermi particles (with two-body interaction) practically does not depend on a 
particular many-body system and has a universal form that can be approximately described 
by the Fermi-Dirac distribution with renormalized parameters (i.e. orbital energies, chemical 
potential, and temperature) 55|-|57|. Statistical effects of the interaction have been shown 
to be absorbed by introduction of the effective temperature. 



In view of the close relationship between strong static correlation effects and representa- 
tions of the ground-state density as well as the close relationship between the distribution 
of occupation numbers for interacting electrons (in the sense of statistical average) and the 
Fermi-Dirac distribution for noninteracting electrons, in this work, the ground-state den- 
sity p(r) of a system of N interacting electrons (at zero temperature) in the presence of 
an external potential v ex t(r), is assumed to be noninteracting thermal ensemble (NI-TE) 
f s -representable, as it is represented by the thermal equilibrium density of an auxiliary 
system of N noninteracting electrons at a fictitious temperature 9 = kBT e i (where ks is 
the Boltzmann constant, T e ; is the temperature measured in absolute temperature, and 9 
is the temperature measured in energy units) in the presence of a local potential v s (r). 
Correspondingly, p(r) can be expressed as 

oo 

p(r) = $>hMr)| 2 , (7) 

i=i 

where the occupation number is the Fermi-Dirac function 

f l = {l+ex P [(e l - f i)/9]}- 1 , (8) 
which obeys the following two conditions, 

oo 

52fi = N, 0</i<l, (9) 

£j is the orbital energy of the i th orbital ipi{v), and p is the chemical potential chosen to 
conserve the number of electrons N. 

In the following section, we demonstrate how a DFT associated with the NI-TE v s - 
representable p(r) in Eq. (JTj), can be formulated. In other words, for a given fictitious 
temperature 9, the remaining " renormalized parameters" ({e^} and p) and the {^(r)}, can 
be self-consistently determined to represent the p(r) in Eq. (J7|), which then determines the 
ground-state energy of the system. Strong static correlation is shown to be described by a 
term related to the 9 and {f{\ in this DFT. 

To avoid any possible confusion with KS-DFT, E-DFT, and finite-temperature DFT 



(FT-DFT) 58], we refer to this DFT as thermally-assisted-occupation DFT (TAO-DFT). 



We wish to develop TAO-DFT with the following characteristics. 
• It is developed for ground-state systems at zero temperature. 



It represents the ground-state density with orbitals and occupation numbers. 

It may be used together with existing XC functionals in KS-DFT. 

It reduces to KS-DFT in the absence of strong static correlation effects. 

It treats single- and multi-reference systems in a more balanced way than KS-DFT. 

It has similar computational cost as KS-DFT (e.g. energy, analytical nuclear gradi- 
ents). 

III. TAO-DFT 

A. Self-Consistent Equations 

Consider a system of N interacting electrons moving in an external potential v ex t(r) at 
zero temperature. Based on the HK theorems the ground-state energy E[p], a functional 
of the ground-state density p(r), can be written as 

E[p] = F[p] + J p(r)v ext (r)dr, (10) 

where the universal functional 

F[p\=T[p\ + V ee [pl (11) 
is the sum of the interacting kinetic energy T[p] and the electron-electron repulsion energy 

In KS-DFT l|, |2j, F[p] is usefully partitioned as 

F[p] = T s [p] + E H [p\ + (T[p] + V M \p\ - T s [p] - E H [p\) 

= T s [p]+E H [p} + E xc [p]. (12) 

Here T s [p] is the noninteracting kinetic energy at zero temperature, 

EM ^jJJ ^« r < (13) 

is the Hartree energy, and E xc [p] = (T[p\ + V ee [p] — T s [p] — E H [p}) is the XC energy defined 
in KS-DFT. In KS-DFT, T s [p], the big unknown in terms of the density, is exactly treated 
by the use of KS orbitals. However, as previously discussed, the basic ansatz of KS-DFT 



(i.e. NI-PS ivrepresentability of the given p(r)) can be violated for SC systems, in which 
case even the exact KS-DFT will fail to convey reliably accurate results [7|. 

To make progress, a different representation of the ground-state density is adopted in 
TAO-DFT, wherein p(r) is represented by the thermal equilibrium density of an auxiliary 
system of N noninteracting electrons at a fictitious temperature 9 in the presence of some 
local potential v s (r). Aiming to achieve this representation, in contrast to the original KS 
partition, F[p] is partitioned into the following set of terms: 



F[p] = A'\p] + E H [p] + (T[p\ + V ee [p] - A» s [p\ - E H [ 

= A 6 s [p] + E H [p] + [T\p] + V ee [p] - T s [p] - E H [p\) + [T s [p] - A e s [p]) 

= A 6 s [p] + E H [p] + E xc [p] + E e [p]. (14) 

Here T s [p], Eh[p], and E xc [p] are the same as those defined in KS-DFT, A e s [p] is the nonin- 
teracting kinetic free energy at temperature 9, and Eq[p] = T s [p] — A s [p] = A e s =0 [p] — A e s [p] is 
the difference between the noninteracting kinetic free energy at zero temperature and that 
at temperature 9. 

Substituting Eq. (Ti~4")) into Eq. (fTUj) and minimizing the E[p] with respect to p(r) (subject 
to the constraint that the number of electrons be N), yields the following Euler equation 
for the ground-state density p(r), 

„ = + „„ ((r) + ? I <^L d , + 'KM + *W (15) 

op(r) J |r — r| op( r ) op(r) 

where p is the chemical potential of the system. 

To bypass the exact functional form of A e s [p] (the big unknown in terms of the density) 
needed in Eq. (1151) , consider an auxiliary system of noninteracting electrons moving in a 



local potential v s (r) at temperature 9. Based on Mermin's theorems 58], the grand-canonical 
potential Q e s of this reference system, a functional of the thermal equilibrium density p s (r), 
can be written as 

tfM = K\Ps] + f Ps(r)[v s (v) - p s ]rfr, (16) 

where p s is the chemical potential of the reference system. 

Minimization of the Q e s [p s ] with respect to p s (r), gives the following Euler equation for 
the thermal equilibrium density p s (r), 

»s = 1 ^y + v s (r). (17) 
10 



Comparing Eq. ( ITT)) with Eq. (TT5)) . shows that both minimizations have the same solution 
Ps(t)=p(t), if we choose v s (r) (up to a constant) as 

Alternatively, as Af[p] can be expressed exactly in terms of orbitals and occupation 
numbers (see below), Eq. (j!7p can also be handled, in an exact manner, by solving the 
one-electron Schrodinger equations for the potential v s (r), given by 

{-^-V 2 + v s (r)}^(r) = e^(r), (19) 



<e 



and construct 

oo 

p(r)= Ps (r) = J2m(r)\ 2 , (20) 

i=i 

where the occupation number /j is the Fermi-Dirac function 

f l = {l + ex P [(e l -p)/9]}- 1 , (21) 
and the chemical potential p is chosen to conserve the number of electrons N, 

oo 

^{l+eMiei-^/Oir^N. (22) 

The formulation of TAO-DFT leads to a set of self-consistent equations, Eqs. ( TT8j) . ( TT9)) . 
(EUD, fl2JJ, and ((22]). 

To obtain a self-consistent ground-state density in TAO-DFT: (i) Choose a trial p(r) to 
construct f s (r) by Eq. (TT8jk (ii) solve Eq. (TT9)) . which gives {ej,^j(r)}; (iii) find /x by solving 
Eq. ( 122)) ; (iv) determine by Eq. (121)) and new p(r) by Eq. ( 120)) . This process is coupled 
with Eq. (118)) to achieve self-consistency. When converged, the entropy functional reads 

oo 

S 6 s [{fi}} = -k B $> ln(/ 4 ) + (1 - /*) ln(l - /«)]. (23) 

i=l 

The exact noninteracting kinetic free energy A s at the fictitious temperature 6, can be 
expressed in terms of {fi,ipi}: 

A e s [{fM] = T^U,^}] - l~ B S 6 s [{h}i (24) 
11 



which is the sum of the kinetic energy 



h 2 



e i=i 



i=l 



p(r)v s (r)dr 



and entropy contribution 



(25) 



i=i 



(26) 



of noninteracting electrons at temperature 9. Based on Eqs. (110]) and (I14j) . the total ground- 
state energy can be evaluated by 

£[p] = A°[{fi, ^}] + y p(r)^(r)dr + £ H [p] + ^ c [p] + W (27) 

To sum up, in TAO-DFT, the partition of F[p] in Eq. (JT3J) and the exact treatment 
of A s [p] in Eq. (|2"%|) . are shown to yield a set of self-consistent equations for the NI-TE v s - 
representable p(r) in Eq. (120]) . Note that these equations resemble the finite-temperature KS 
equations 2j, [58|, so the implementation of TAO-DFT can be easily achieved using existing 
FT-DFT codes with a slight modification (i.e. replacing the XC free energy in FT-DFT to the 
sum of E xc [p] and E$[p] in TAO-DFT). Hence, the computational cost of TAO-DFT is similar 
to that of KS-DFT or FT-DFT. Similar to FT-DFT [2, 5f|, due to the explicit inclusion of 
Fermi-Dirac occupation function in TAO-DFT, the entropy contribution (— -^-S e s \{fi\\] 
Eq. (126]) is essential to make the total ground-state energy functional E[p] variational 
(e.g. making the nuclear gradients of E[p] equal to the Hellmann-Feynman forces jfiQl). 
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B. Spin-Polarized Formalism 

For a system with N a up-spin electrons and Np down-spin electrons, the standard compu- 
tational approach is the spin-polarized (spin-unrestricted) formalism, wherein the fundamen- 
tal variables are the up-spin density p a (r) and down-spin density pp(r) of the ground-state 
density 

p{r) = p a (r) + p p {v) = £ p CT (r). (28) 

a=a,fi 

12 



In analogy to the two- Fermi- level picture of spin-polarized KS-DFT 48|, |6l), spin-polarized 
TAO-DFT can also be formulated with the two-chemical-potential picture, wherein two non- 
interacting auxiliary systems at the same fictitious temperature 9 are adopted: one described 
by the spin function a and the other by function 0, with the corresponding thermal equilib- 
rium density distributions p s , a (r) and p s ,p{r) exactly equal to p a (r) and pp(r), respectively, 
in the original interacting system at zero temperature. Similar to the previous derivations 



(but using the spin-polarized extensions of the HK theorems 48|, |6lfl and the Mermin theo 



rems 



62 



63| for the physical and auxiliary systems, respectively), one-electron Schrodinger 
equations for electrons with cr-spin (cr = a or (3), can be obtained as follows (i runs for the 
orbital index): 



h 2 



2m, 



-V 2 + u», ff ( r )}Vv( r ) = e i,<xVv( r )> 



with a local potential 



Jext\ 



+ e 2 f P(r') dr , + 5E xc [p a , Pl3 } + SE e [p a , 
J |r - r'| 8p<r(r) Sp a (i 



(29) 



(30) 



Here E xc [p a , pp] is the same as the XC energy defined in spin-polarized KS-DFT 48|, |6jJ], and 



EeiPaiPp) = T„\p ai pp] - A e s [p a ,pp\ = A e s °{p a ,pf3] - A e 8 \p a ,p p ] is the spin-polarized version 
of Eq[p]. The cr-spin density can be constructed by 



i=i 

where the occupation number fa a is the Fermi-Dirac function 



(31) 



fi,a = {1 + exp[(e ij(T - Pa)/0]} 



-1 



(32) 



and a chemical potential p a is chosen to conserve the number of cr-spin electrons N a , 

oo 

^{1 + exp[(e li(T - p^/9]}- 1 = N a . (33) 



The formulation of spin-polarized TAO-DFT yields two sets (one for each spin function) of 
self-consistent equations, Eqs. (129 p . (130 p . ( I3ip . ( I32p . and (I33p . for p a (r) and pp(r), respec- 
tively, which are coupled with p(r) by Eq. (|28j) ■ 

To obtain self-consistent spin densities (and the ground-state density) in spin-polarized 
TAO-DFT: (i) Choose trial spin densities p a (r) and pp{v) to compute the ground-state 
density p(r) by Eq. (1281) : (ii) for the cr-spin (a = a or 0) electrons, construct v s>a (r) by 



13 



Eq. ([H; (iii) solve Eq. ([29]), which gives {e ii<T , ^(r)}; (iv) find p ff by solving Eq. ([33]); (v) 
determine {/i, CT } by Eq. ([32]) and new Po-(r) by Eq. ([3T]) . This process is coupled with Eq. 
([28]) to achieve self-consistency. When converged, A e s a , the sum of the kinetic energy and 
entropy contribution of noninteracting a-spin electrons at the fictitious temperature 9, is 
given by 

fc2 00 /* 00 

< CT [{ A., W}] = E A* / ^(r) V 2 ^,.(r)rfr + In(/^) + (1 - m(l - / ii<7 )] 

e i=l " i=l 

oo „ 

= ^{fi^ + [/ i>ff ln(/ ijCT ) + (1 - fi t a) ln(l - - / p ff (r)tv(r)dr, (34) 

i=i ^ 

and the total ground-state energy ^[p^p^] in spin-polarized TAO-DFT is evaluated by 

#[Pa,P/3] = E ^1<t[{./W>Vv}]+ / P( r ) V e*t(r)tfr+^H[p] + #*c[Pa, P/j] + ^[Pa, P/j]. (35) 
<r=a,/3 

Spin-unpolarized (spin-restricted) TAO-DFT can be formulated by imposing the con- 
straints of xj)i j0l (r) = ip it p(r) and fi a = f^p to spin-polarized (spin-unrestricted) TAO-DFT. 



C. Analytical Nuclear Gradients 

The analytical computation of nuclear gradients is crucial for the efficient optimization of 
molecular geometries. In light of the similarity of TAO-DFT and FT-DFT 0,0], analytical 
nuclear gradients for TAO-DFT can be easily obtained from those for FT-DFT with a slight 
modification (mentioned previously). Therefore, the computational cost of the analytical nu- 
clear gradients for TAO-DFT is similar to that for KS-DFT or FT-DFT. For nonorthogonal 
atomic-orbital representations (e.g. Gaussian-type orbitals), the generalized Pulay force for 



a noninteracting thermal ensemble 



64] should be included to add the effect of the basis-set 



dependent response to the analytical nuclear gradients for TAO-DFT. 



D. Local Density Approximation 

In TAO-DFT, as the exact E xc [p] and Eg[p], in terms of the ground-state density p(r), 
remain unknown, DFAs for both of them (denoted as TAO-DFAs) are needed for practical 
applications. The performance of TAO-DFAs depends on the accuracy of DFAs and the 
choice of the fictitious temperature 9. In this work, we adopt the LDA (the simplest DFA) 
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for both the E xc [p] and E e [p] in TAO-DFT (denoted as TAO-LDA). As TAO-LDA is exact 
for a uniform electron gas (UEG), it provides a good starting point for more accurate and 
sophisticated TAO-DFAs. Besides, TAO-LDA is readily available, as E x jBA [p] can be directly 



obtained from that of KS-LDA 
of A^ DA ' e [p\ as follows: 



65 



66fl, and E% DA [p\ can be obtained with the knowledge 



E^ A [p]=A 



LDA,0=O 



[p]-A^[p]. 



(36) 



Here, Perrot's parametrization of ' [p\ (in its spin-unpolarized form) 67J] is adopted to 



by Eq. f)36p . For completeness of this work, 



obtain Eq [p] (in its spin-unpolarized foruf 

v4^ DA ' e [p] (after correcting some typos in Ref. 67]) is explicitly shown here (in atomic units) 



.4 



LDA,6» 



[p] 



,LDA,i 



(r)dr, (37) 



unction f(y) was parametrized 



67J: 



where a^ DA ' e (r) = 6p(r)f(y) andy = (n 2 /^2)0- 3 / 2 p(r). The: 
separately for the two regions y < y and y > y (y = j^jg) 

f(y) = \ny - 0.8791880215 + 0.1989718742?/ + 0.1068697043 x 10" V 

- 0.8812685726 x I0~ 2 y 3 + 0.1272183027 x lO^y 4 - 0.9772758583 x 10~V 

+ 0.3820630477 x lO" 2 ?/ 6 - 0.5971217041 x 10~V, for V < Vo, (38) 



f(y) = 0.7862224183m - 0.1882979454 x lO 1 ^ 1 + 0.5321952681m"' 3 

+ 0.2304457955 x lO 1 ?^ 5 - 0.1614280772 x 10 2 u" 7 + 0.5228431386 x 10V 

- 0.9592645619 x 10 2 M~ n + 0.9462230172 x 10 2 m~ 13 

- 0.3893753937 x 10V 15 , for y > y with u = y 2/? \ 



(39) 



The 9 = case, A^ BA,e °[p\, is the same as the Thomas- Fermi kinetic energy density func- 
tional {3, 4|, 



.4 



LDA,6»=0r 



[p]=C F I p 5 /\v)dv 



(40) 



3 ; (3vr 2 ) 2 / 3 . 



where Cp — < 

For spin-polarized (spin-unrestricted) TAO-LDA, the corresponding spin-polarized forms, 



65 



66j) and E^ BA [p a 



E%c [p a , Pp] (available from that of spin-polarized KS-LDA 
should be adopted. From the spin-scaling relation of A e s [p a , pp] (same as that of T s [p a ,pp] 
68]), E^ DA [p a , pp] can be conveniently expressed by its spin-unpolarized form E^ BA [p], 

E^ A [ Pa , P p] ee A^ s =°[p a ,pp] - A^ A ' e [ Pa ,pp] 
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\{A™^[2 Pa \ + A^°[2p,\) - \{A^[2 Pa ] + A^[2p,]) 
\{{A^^[2 Pa ] - A^[2p a }) + (A^=°[2 Pp ) - A^[2 P ,])} 
l -{E^ A [2p a } + E^[2p,]}. (41) 



E. Strong Static Correlation from TAO-LDA 



The ground-state density p(r) of a strongly correlated system containing a sufficiently 
large number of electrons, can be represented by Eq. §T§ (with the exact NOs {Xi{ r )} an d 
NOONs {rii}). Assume that the p(r) can also be represented by Eq. (120]) (with the orbitals 
{"0i(r)} and their occupation numbers {fi} from the exact TAO-DFT). Note that for such a 
NI-TE w s -representable p(r), its "internal variables" {fi} and {^(r)} in Eq. (T2T)j) . can still 
be tuned by changing the fictitious temperature 9. If a 9 is chosen so that {n{} ~ {fi} (in 
the sense of statistical average, as mentioned previously), we have {Xi{ r )} ~ {^M 1 ")} ( as 
both Eq. ([1]) and Eq. (1201) represent the same p(r)). In fact, the NI-TE t> s -representability 
of p(r) is likely to be fulfilled for this 9, due to the similarity of Eq. (pQ) and Eq. (120|) . 

Consequently, the exact kinetic energy of interacting electrons T[p] can be properly sim- 
ulated by Tf [{f,i>i}] (as appeared in A a [{fi,if>i}] = T e s [{f i} ipi}] - ^S e s [{fi}]), namely, 

T[p\=T[{n uXl }}^T e s [{f l ^ l }} 1 (42) 

due to their similar expressions j^j, while the electron-electron repulsion energy V ee [p] = 
F[p] - T[p\ (see Eqs. (PU and (HH)) is given by, 

V ee [p] « E H [p] + E xc [p] + Eeip] - ^S e s [{fi}}. (43) 

On the right-hand side of Eq. (l4~3]) . the first term is the Hartree energy, and the sum of 
the remaining terms (E xc [p\ + Eg[p] — -^-Sf [{fi})) should properly describe the XC energy 
defined in the exact wave function theory. 

Here, we explain how strong static correlation is described by TAO-LDA, with arguments 
similar to the above. Suppose that the exact p(r) in Eq. (CD) can be reasonably represented 
by Eq. ( 1201) (with the orbitals {^(r)} and their occupation numbers {f} from TAO-LDA). 
When applying the above arguments for TAO-LDA (i.e. choosing a 9 so that {n^} ~ {fi}, 
which gives {xi(r)} ~ {ipi(r)}), T[p] can still be properly simulated by T s e [{/j, ipi}] (see Eq. 
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2J)), while V ee [p] is only approximated by, 



V M \p] « ShW + ^ L C DA [P] + ^ LDA [P] - ^StM}). (44) 

On the right-hand side of Eq. (T4*4"j) . the first two terms (-E/f [p] and E^ A [p\) are the same as 
those defined in KS-LDA, the third term E^ DA [p] locally accounts for the difference between 
the exact T s and A e s (at the LDA level), and the last term is the entropy contribution 
(-&#[{/<}] « -£S*[{ni}] = ^E^iN ln(n,) + (1 - n,) In(l - m)]) (see Eq. ». 

Due to their local treatment, i?^? [p] and i?g DA [p] are not expected to properly de- 
scribe nonlocal XC effects (e.g. long-range dynamical correlation, strong static correlation). 
However, as the entropy contribution is a fully nonlocal density functional ({fi} are implicit 
density functionals), it may describe nonlocal correlation effects. There is certainly a close re- 
lationship between the entropy (defined by the NOONs {ni}) and correlation energy of a sys- 



tem. A famous example is given by the Collins conjecture 



691 ] that the correlation energy of 



a system is proportional to the Jaynes (information) entropy Sj a ynes[{^i}] — — n i ^ n ( n i) 



70j. Interestingly, the entropy contribution in Eq. (|44p is proportional to the Gibbs (ther- 



modynamic) entropy (Sf[{/i}] « ^[{ni}] = -A;bZ)£iK m (^i) + (1 - th) ln(l — ra*)]), with 
the constant of proportionality being explicitly given by (— -j^). Note that the similarity 
of information entropy and thermodynamic entropy in a many-body quantum system (with 
strong interactions) has been shown in Ref. 7l|, based on statistical arguments. 

As the entropy contribution (— -^S d s [{fi})) in Eq. (T4*4|) essentially provides no contribu- 
tions for a single-reference system ({fi} ~ {rii} are close to either or 1), and significantly 
lowers the total energy of a multi-reference system ({fi} ~ {n,i\ are fractional for active 
orbitals, and are close to either or 1 for others), we expect that this term (absent in 
KS-LDA) play a crucial role in simulating strong static correlation (rather than dynamical 
correlation) . 



IV. NUMERICAL INVESTIGATIONS OF AN OPTIMAL 9 VALUE 

The fictitious (reference) temperature 9 for TAO-DFT, controlling the orbital occupation 
numbers {fi}, is closely related to the strength of static correlation. At the LDA level, an 
immediate question is how the 9 for the resulting TAO-LDA should be chosen. As previously 
argued, the entropy contribution, -^Sf [{/»}] = 9 YnLiifi ln (/i) + (1 ~ fi) ln (l ~ fi)), can 
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be responsible for strong static correlation effects, especially when the {f{\ (tunable by the 
9) properly simulate the exact NOONs {rii}. To numerically investigate this conjecture, 
the performance of TAO-LDA (with 9 = 0, 1, 3, 5, 7, 10, 15, 20, 30, and 40 mHartree) is 
examined for both single-reference systems (reaction energies and equilibrium geometries) 
and multi- reference systems (dissociation of H 2 and N 2 , twisted ethylene, and singlet-triplet 
energy gaps of linear acenes). The limiting case where 9 = for TAO-LDA is especially 
interesting, as this reduces to KS-LDA. Therefore, it is important to know how well KS-LDA 
performs here to assess the significance of the extra parameter 9 for TAO-LDA. 



All calculations are performed with a development version of Q-Chem 3.2 72]. The error 
for each entry is defined as (error = theoretical value — reference value). The notation used 
for characterizing statistical errors is as follows: mean signed errors (MSEs), mean absolute 
errors (MAEs), root-mean-square (rms) errors, maximum negative errors (Max(— )), and 
maximum positive errors (Max(+)). Results are computed using the 6-311++G(3df,3pd) 
basis set, unless noted otherwise. 



A. Single-Reference Systems 



1. Reaction Energies 



The accurate prediction of reaction energies is usually one of the major criteria in the 
assessment of the performance of electronic structure methods. The reaction energies of 30 



chemical reactions (a test set described in Ref. |20j) are used to examine the performance of 
TAO-LDA. As shown in Tabled! TAO-LDA (with a 9 smaller than 10 mHartree) has similar 
performance to KS-LDA 7jJ. This is unsurprising, as these systems do not have much static 
correlation, the exact NOONs should be close to either or 1, which can be well simulated 
by the orbital occupation numbers of TAO-LDA (with a sufficiently small 9). Consequently, 
T ![{fi^i}} ( see Eq. (E2D) is close to Tf =0 [{^}] (KS kinetic energy), and E^ DA [p] (see Eq. 
(J36D) and the entropy contribution (-££*[{/,}] = 9J2T=Afi Hfi) + i 1 ~ fi) ln (l - /<)]) 
have insignificant contributions to the total energy, relative to E^ A [p] (see Eq 
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2. Equilibrium Geometries 



Geometry optimizations for TAO-LDA are performed on the equilibrium experimental 



test set (EXTS) 



74J ]. consisting of 166 symmetry unique experimental bond lengths for small 



to medium size molecules. As the ground states of these molecules near their equilibrium 
geometries can be well described by single-reference wave functions, TAO-LDA (with a 9 
smaller than 10 mHartree) is also found to perform similarly to KS-LDA {7^, as shown in 
Table HI 



B. Multi- Reference Systems 

1. Dissociation of Hi and N2 

H 2 dissociation, a single-bond breaking system, is particularly challenging for KS-DFT. 
Fig. (JT]) shows the potential energy curves (in total energy) for the ground state of H 2 , cal- 
culated by both the s pin -restricted and spin-unrestricted formalisms of the HF theory and 
KS-DFT (with LDA jfifi, fifj 



and B3LYP 



13 



14j functionals), where the exact potential 



energy curve is calculated by the CCSD theory (coupled-cluster theory with iterative singles 
and doubles) 75|. Due to the symmetry constraint, the spin-restricted and spin-unrestricted 
potential energy curves, calculated by the exact theory, should be the same. Therefore, the 
difference between the dissociation limits of the spin-restricted and spin-unrestricted poten- 



tia 



10 



energy curves, can be used as a quantitative measure of SCEs of approximate methods 
, |35j. Spin- restricted KS-DFT yields the proper spin symm etry and spin densities, but 
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due to the lack of 



has much too high total energy (leading to a noticeable SCE 
strong static correlation. On the other hand, spin-unrestricted KS-DFT artificially breaks 
the correct space- and spin-symmetries to simulate strong static correlation, yielding a rea- 
sonable energy but wrong spin densities 10j. Similar results are also found for the HF theory. 



Among the three approximate methods, HF has the largest SCE due to the complete neglect 
of static (and also dynamical) correlation. The SCE of LDA is smaller than that of B3LYP 
or HF. In Fig. (j2J), the potential energy curves (in relative energy) for the ground state 



of H2, calculated by the spin-restricted HF theory and KS-DFT 



3 



(with LDA 



76 , B3LYP 



13 



14|], M06-2X [15], wB97X-D |24j, and B2PLYP 



65j,l66j, PBE 



3l| functionals), are pre- 



sented for comparisons, where the zeros of energy are set at the respective spin-unrestricted 
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dissociation limits. As can be seen, LDA still has the smallest SCE, when compared with 
other approximate methods. Popular hybrid functionals (B3LYP, M06-2X, and cuB97X-D) 
perform very well near the equilibrium geometry (dominated by single-reference character), 
but fail drastically at the larger R (dominated by multi- reference character). B2PLYP (a 
popular DH functional) leads to an unphysical divergence at the dissociation limit, due to 
the vanishing HOMO-LUMO gap appeared in the energy denominator of its second-order 
perturbation energy components. Clearly, hybrid and DH functionals, the most popular 
schemes for reducing the SIEs and NCIEs of KS-DFAs, respectively, can perform poorly for 
multi-reference systems due to their inaccurate treatment of strong static correlation effects 



10 



35j|. 



To evaluate the performance of the present method, the potential energy curves for the 
ground state of H2, calculated by spin- restricted TAO-LDA, are shown in Fig. [3] (in total 
energy) and Fig. H] (in relative energy). Near the equilibrium geometry, where the multi- 
reference character is insignificant, the performance of TAO-LDA (with a 9 smaller than 10 
mHartree) is very similar to that of KS-LDA. At the dissociation limit, where the multi- 
reference character is pronounced, the SCE of TAO-LDA is shown to be reducible with the 
increase of 9 value, at essentially no extra computational cost! Interestingly, TAO-LDA 
(with a 9 between 30 and 50 mHartree) performs very well, leading to a vanishingly small 
SCE. 

To see how this is related to the ensemble representation (via fractional orbital occu- 
pations) of the ground-state density, the occupation numbers of the \a g orbital (HOMO) 
for the ground state of H2 as a function of the internuclear distance R, calculated by spin- 
restricted TAO-LDA, are presented in Fig. El where the reference data are the FCI NOONs 
0. At the equilibrium geometry (R = 0.741 A), the FCI NOON is 1.9643, indicating the 
absence of strong static correlation effects (with doubly occupied \a g orbital) . However, the 
FCI NOON is 1.5162 at R = 2.117 A, and 1.0000 at R = 7.938 A, indicating the presence 
of strong static correlation effects. The la g orbital occupation numbers of spin-restricted 
TAO-LDA (with a nonvanishing 9) are very close to 2.0 (doubly occupied) near the equilib- 
rium geometry, and gradually reduced to 1.0 (singly occupied) at the dissociation limit. The 
larger the 9 value is, the faster the corresponding la g orbital occupation number approaches 
1.0 at the larger R. The \a g orbital occupation numbers of spin-restricted TAO-LDA (with 
a 9 between 30 and 50 mHartree) are shown to match well with the FCI NOONs, which is 
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closely related to the vanishingly small SCE of TAO-LDA (with the same 9). 

To examine the entropy contributions (in total energy) as a function of the internuclear 
distance R, calculated by spin- restricted TAO-LDA (6 = 40 mHartree), Fig. shows the 
potential energy curves (in relative energy) for the ground state of H2, calculated by the spin- 
restricted (with and without the entropy contributions) and spin-unrestricted TAO-LDA (9 
= 40 mHartree), where the zeros of energy are set at the spin-unrestricted dissociation limit. 
As can be seen, the entropy contributions are insignificant near the equilibrium geometry 
and pronounced at the larger R (approaching a negative constant at the dissociation limit), 
properly simulating the strong static correlation effects to make the spin-restricted potential 
energy curve the same as the spin-unrestricted one (as it should be). By contrast, the entropy 
contributions as a function of the internuclear distance R, calculated by spin-restricted TAO- 
LDA (9 = 7 mHartree), are still insufficient to simulate the strong static correlation effects 
(see Fig. [7]), as the corresponding la g orbital occupation numbers do not match well with 
the FCI NOONs (see Fig. EJ. 

Similar results are also found for N 2 dissociation, a triple-bond breaking system. The 
potential energy curves for the ground state of N 2 , calculated by spin- restricted TAO-LDA, 
are shown in Fig. [8] (in total energy) and Fig. [9] (in relative energy). As can be seen, spin- 
restricted TAO-LDA (with a 9 between 30 and 50 mHartree) can dissociate N 2 properly 
(yielding a vanishingly small SCE) to the respective spin-unrestricted dissociation limits, 
which is closely related to that the occupation numbers of the 3a g (in Fig. [TO]) and lft ux 
(in Fig. [TT]) orbitals for the ground state of N 2 as functions of the internuclear distance R, 
calculated by spin-restricted TAO-LDA (with the same 9), match reasonably well with the 
corresponding MRCI NOONs (the reference data) jz7|. 

To examine the entropy contributions (in total energy) as a function of the internuclear 
distance R, calculated by spin-restricted TAO-LDA (9 = 40 mHartree), Fig. IT2l shows the 
potential energy curves (in relative energy) for the ground state of N 2 , calculated by the spin- 
restricted (with and without the entropy contributions) and spin-unrestricted TAO-LDA (9 
= 40 mHartree), where the zeros of energy are set at the spin-unrestricted dissociation limit. 
As shown, the entropy contributions are essentially responsible for simulating the strong 
static correlation effects to make the spin-restricted potential energy curve the same as the 
spin-unrestricted one (as it should be). For spin-restricted TAO-LDA (9 = 40 mHartree), 
the entropy contribution (-208.90 kcal/mol) at the dissociation limit of N 2 , is considerably 
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larger (about three times larger) than that (-69.64 kcal/mol) at the dissociation limit of H2, 
as the number of unpaired electrons (or singly occupied orbitals) for N2 dissociation is more 
than (three times more) that for H 2 dissociation. 

To sum up, when the orbital occupation numbers {/j} of TAO-LDA are close to the exact 
NOONs {rii}, the strong static correlation effects are shown to be properly simulated by the 
entropy contribution of TAO-LDA. As this feature is independent of the number of unpaired 
electrons in a system, TAO-LDA seems to be promising for the study of large polyradical 
systems, such as linear acenes (as will be shown later). 



2. Twisted Ethylene 

The 7r (lb 2 ) and ir* (2b 2 ) orbitals in ethylene (C2H4) should be degenerate when the 
HCCH torsion angle is 90°. Spin-restricted single- reference methods cannot handle such a 
degeneracy properly and show an unphysical cusp in the torsion potential near 90°. In the 
calculations, we use the experimental geometry of C2H4 (Rcc = 1-339 A, Rqu = 1.086 A, 
^hch = 117.6°) [78|. Fig. [TBI shows the torsion potential energy curves (in relative energy) 
for the ground state of twisted ethylene as a function of the HCCH torsion angle, calculated 
by spin-restricted TAO-LDA, where the zeros of energy are set at the respective minimum 
energies. Spin-restricted TAO-LDA (with a 9 larger than 5 mHartree) is shown to be able 
to remove the unphysical cusp, though spin-restricted TAO-LDA (with a 9 larger than 20 
mHartree) is shown to yield a torsion barrier which is far too low. 

Fig. [T41 shows the occupation numbers of the 7r (lb 2 ) orbital for the ground state of twisted 
ethylene as a function of the HCCH torsion angle, calculated by spin-restricted TAO-LDA, 
where the reference data are the half-projected NOONs of CASSCF method (HPNO-CAS) 



79|. As can be seen, the ir (lb 2 ) orbital occupation numbers of spin-restricted TAO-LDA 
(with a 9 between 10 and 20 mHartree), match reasonably well with the accurate NOONs, 
which is related to the accurate torsion potential energy curve, calculated by spin-restricted 
TAO-LDA (with the same 9). 

To examine the entropy contributions (in total energy) as a function of the HCCH torsion 
angle, calculated by spin- restricted TAO-LDA {9 = 15 mHartree), Fig. H~5l shows the torsion 
potential energy curves (in relative energy) for the ground state of twisted ethylene as a 
function of the HCCH torsion angle, calculated by the spin-restricted (with and without the 
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entropy contributions) and spin-unrestricted TAO-LDA (8 = 15 mHartree), where the zeros 
of energy are set at the respective minimum energies. As shown, the entropy contributions 
are responsible for simulating the strong static correlation effects to make the spin-restricted 
torsion potential energy curve the same as the spin-unrestricted one (as it should be). 



3. Singlet-Triplet Energy Gaps of Linear Acenes 



Linear n-acenes (C4 n+ 2H 2n +4), consisting of n linearly fused benzene rings (see Fig. (fTfijP 



have attracted great interest from many experimental and theoretica 
their fascinating electronic properties and technological potential 



researchers due to 



80|-|95|. The experimental 



80 



830, 



singlet-triplet energy gaps (ST gaps) of n-acenes are only available up to pentacene 
due to the increasing reactivity of the larger acenes. Recently, the calculated ST gaps have 
been in serious debate 



84j-l94j. Typically, delocalized 7r-orbital systems, such as n-acenes, 



require high-level ab initio methods, such as the DMRG algorithm 



88] . to capture the 



essential strong static correlation effects. Based on the recent work of Chan and co-workers 



88j, the DMRG ST gaps as a function of the acene length have been shown to decrease 



monotonically with increasing chain length. Based on a good fit to the DMRG ST gaps of 
the smaller n-acenes (up to 12-acene), an exponential fitting function of the form a + b e~ c n 
was adopted for extrapolation of the ST gaps to the infinite chain limit [88], yielding a 
finite ST gap (3.33 kcal/mol for the cc-pVDZ basis set, and 8.69 kcal/mol for the STO-3G 
basis set) for polyacene (triplet above singlet). However, the extrapolated results have been 



shown subject to details of the fit 



92j. More importantly, it is unclear whether the ST 



gaps of the larger n-acenes (eg. n > 20) still decrease exponentially (as those of the smaller 
n-acenes) with increasing chain length, which may significantly affect the extrapolated ST 
gap. Calculations on the larger acenes are necessary to address this, which are, however, 



prohibitively expensive for the DMRG algorithm 
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881 ] and other high-level ab initio methods 



On the other hand, KS-DFT is computationally efficient, but unable to handle such strong 



static correlation effects properly. To show t 



LDA 
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BLYP 



96, 



lis, spin-unrestricted KS-DFT calculations (with 



and B3LYP 



13 



14j functionals) are performed, using the 6- 



31G* basis set (up to 16-acene), for the lowest singlet and triplet energies on the respective 
geometries that were fully optimized at the same level. The ST gap of n-acene is calculated 
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as {Ei — Es), the energy difference between the lowest triplet (T) and singlet (S) states of 



n-acene. As shown in Fig. (|17j) . in contrast to the DMRG results 88], the ST gaps calculated 
by spin-unrestricted KS-DFT, are shown to unexpectedly increase beyond 10-acene, due to 
unphysical symmetry-breaking effects 

To assess the performance of the present method, spin-unrestricted TAO-LDA calcula- 
tions are performed, using both the 6-31 G* (up to 46-acene) and 6-31 G (up to 74-acene) 
basis sets, for the lowest singlet and triplet energies on the respective geometries that were 
fully optimized at the same level In Fig. f[T8"|) . the calculated ST gaps as a function of 
the acene length (using the 6-31G* basis set) are plotted. In contrast to the spin-unrestricted 
KS-LDA results, the ST gaps calculated by spin-unrestricted TAO-LDA (with a 9 larger than 
5 mHartree), are shown to decrease monotonically with the increase of chain length. The ST 
gaps calculated by spin-unrestricted TAO-LDA (with a 9 between 5 and 10 mHartree), are 
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in good agreement with the existing experimental and high-level ab initio data (88| 
Due to the symmetry constraint, the spin- restricted and spin-unrestricted energies for the 
lowest singlet state of n-acene, calculated by the exact theory, should be the same. To ex- 
amine this property, spin-restricted TAO-LDA calculations (using the 6-31G* basis set) are 
also performed for the lowest singlet energies on the respective geometries that were fully 
optimized at the same level. The spin-unrestricted and spin-restricted TAO-LDA (with a 9 
larger than 5 mHartree) calculations are found to essentially yield the same energy value for 
the lowest singlet state of n-acene (i.e. no unphysical symmetry-breaking effects). 

Fig. (fT9|) shows the ST gaps of n-acenes (n > 8) as a function of the acene length, 
calculated by spin-unrestricted TAO-LDA {9 = 7 mHartree), using both the 6-3 1G* and 
6-3 1G basis sets {73 1. The effects of basis sets on the calculated ST gaps are shown to 
be insignificant for the larger acenes. The ground state of n-acene is found to remain a 
singlet as chain length is increased. At the level of spin-unrestricted TAO-LDA {9 = 7 
mHartree) /6-31G, the ST gap of the largest acene studied here (74-acene) is 0.66 kcal/mol. 
In view of the slow convergence of the ST gaps with the increase of chain length, the ST 



gaps of the larger n- acenes {n > 20) are found to fit extremely well to a 



power 
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aw function 
I. As shown 



of the form a + b n~ c , rather than the popular exponential function 
in Table IIIIl nonlinear least-squares fitting of 3 different data sets (20- to 74-acene, 30- to 
74-acene, and 40- to 74-acene) of the ST gaps calculated by spin-unrestricted TAO-LDA {9 
= 7 mHartree) /6-31G, by means of the above power-law fitting function, gives estimates of 
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0.08, 0.04, and 0.03 (kcal/mol), respectively, for the ST gaps of n-acenes in the polymer limit 
(n — > oo). As the extrapolated ST gaps are rather insensitive to the choices of the fitting 
data sets, the ST gaps selected in the fitting data sets should have approached the asymptotic 
(large-n) behavior, decreasing as slowly as about rT l with increasing chain length. In view 
of the minor dependence of the extrapolated results with the choices of the fitting data sets, 
we can only conclude that in the polymer limit, the lowest singlet and triplet states should 
be degenerate within 0.1 kcal/mol (triplet above singlet), which supports the absence of 



Peierls distortions 



98] in this limit, and the closure of the fundamental gap 
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The orbital occupation numbers of TAO-LDA provide information useful in assessing the 
possible polyradical character of n-acenes. Fig. (120]) shows the HOMO occupation numbers 
for the lowest singlet states of n-acenes as a function of the acene length, calculated by 
spin-restricted TAO-LDA/6-31G*, where the reference data are the NOONs of the active- 
space variational 2- RDM method [95]. Here, HOMO is the (N/2) th orbital, and LUMO 
is the (N/2 + l) th orbital, where iV is the number of electrons in n-acene. As can be 
seen, the HOMO occupation numbers of spin-restricted TAO-LDA (with a 9 between 5 
and 15 mHartree), match reasonably well with the NOONs, which may suggest that the 
ST gaps calculated by TAO-LDA (with a 9 between 5 and 15 mHartree) should be reliably 
accurate (due to the appropriate treatment of strong static correlation effects via the entropy 
contribution), providing that these agreements are extendible for the larger acenes. 

Fig. (}2~Tj) shows the active orbital occupation numbers for the lowest singlet states of 
n-acenes as a function of the acene length, calculated by spin-restricted TAO-LDA (9 = 
7 mHartree) /6-31G* [73]. Here, for simplicity, HOMO, HOMO-1, and HOMO-6, are 
denoted as H, H-l, and H-6, respectively, while LUMO, LUMO+1, and LUMO+6, 
are denoted as L, L + 1, and L + 6, respectively. As can be seen, the active orbital occu- 
pation numbers exhibit oscillatory behavior in the approach to unity (singly occupied) with 
increasing chain length. The number of factionally occupied orbitals is shown to increase 
with the increase of chain length, which supports previous finding that large acenes should 



exhibit polyradical character 



89|. 



To sum up, it seems plausible to believe the results obtained by TAO-LDA 



= 7 



mHartree) here, as the calculated S' 
mental and high-level ah initio data 



gaps 
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are in good agreement with the existing experi- 



94] . the calculated HOMO occupation numbers 



match reasonably well with the accurate NOONs, and no unphysical symmetry-breaking 
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effects occur for the lowest singlet states of n-acenes. 

V. DEFINITION OF AN OPTIMAL 6 VALUE 

In our study, TAO-LDA (with some fictitious temperature 9) has been found to perform 
reasonably well for multi-reference systems, when the orbital occupation numbers {/j} are 
close to the exact NOONs {n^}. In such a situation, the strong static correlation effects can 
be properly simulated by the entropy contribution of TAO-LDA. However, for multi-reference 
systems, the optimal 9 for TAO-LDA has been found to be highly system-dependent, ranging 
from 5 to 50 mHartree, due to the different strengths of static correlation and the diversities 
of the {rii}. On the other hand, for single-reference systems (in the absence of strong static 
correlation effects), TAO-LDA (with a 9 smaller than 10 mHartree) has been shown to 
perform similarly to KS-LDA. 

For TAO-LDA, although it is impossible to choose a 9 that is optimal for all the systems 
studied, it is still useful to define one to provide an explicit description of orbital occupations. 
Here, the optimal 9 value is defined as the largest 9 value for which the performance of the 
TAO-LDA (with this 9) and KS-LDA is similar for single-reference systems. Based on our 
numerical investigations, an optimal value of 9 = 7 mHartree, is finally chosen. TAO-LDA 
(9 = 7 mHartree) has been shown to consistently improve upon KS-LDA for multi-reference 
systems, while performing similarly to KS-LDA for single-reference systems. 

VI. CONCLUSIONS 

We have proposed TAO-DFT, a DFT with fractional orbital occupations produced by 
the Fermi-Dirac distribution (in order to simulate the distribution of orbital occupation 
numbers for interacting electrons). TAO-DFT offers an explicit description of strong static 
correlation via the entropy contribution, a function of the fictitious temperature 9 and orbital 
occupation numbers {fi} (implicit density functionals). Even at the simplest LDA level, the 
resulting TAO-LDA has been shown to perform reasonably well for multi-reference systems 
(due to the appropriate treatment of static correlation), when the {fi} (related to the 9) are 
close to the exact NOONs {n^}. As this feature is independent of the number of unpaired 
electrons in a system, TAO-LDA seems to be very useful for the study of large polyradical 
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systems. In our study, an optimal value of 9 = 7 mHartree, has been defined based on 
physical arguments and numerical investigations. TAO-LDA (9 = 7 mHartree), though not 
optimal for all the systems studied, has been shown to consistently improve upon KS-LDA for 
multi-reference systems, while performing similarly to KS-LDA for single-reference systems. 
Due to its computational efficiency and reasonable accuracy, TAO-LDA has been applied 
to the study of the ST gaps of acenes, which are challenging problems for conventional 
electronic structure methods. At the level of TAO-LDA (9 = 7 mHartree) /6-31G, the ST 
gap of polyacene has been shown to be vanishingly small (within 0.1 kcal/mol), and large 
acenes should exhibit singlet polyradical character in their ground states. 

As TAO-LDA is conceptually simple, computationally efficient, and easy to implement, 
it seems to be a promising method for the study of ground states of large single- and multi- 
reference systems. However, as with all approximate electronic structure methods, some 
limitations remain. The optimal 9 = 7 mHartree for TAO-LDA is system-independent 
(not fully optimized for each system), a system-dependent 9 (related to the distributions of 
NOONs) is expected to enhance the performance of TAO-LDA for a wide range of systems. 
For single-reference systems, the performance of TAO-LDA is similar to that of KS-LDA. 
A possible TAO-DFA is expected to perform better than TAO-LDA for single-reference 
systems. Although the SCEs of TAO-DFAs are expected to be less than those of KS-DFAs, 
the SIEs and NCIEs of TAO-DFAs may remain enormous in situations where these failures 
occur. A fully nonlocal TAO-DFT (i.e. nonlocal E xc [p] and E e [p\) may be needed to resolve 
all the three qualitative errors (SIE, NCIE, and SCE). We are currently investigating along 
these lines, and results may be reported elsewhere. 
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FIGURES 



FIG. 1. Potential energy curves (in total energy) for the ground state of H2, calculated by both 
the spin-restricted and spin-unrestricted formalisms of the HF theory and KS-DFT (with LDA and 
B3LYP functionals). The exact potential energy curve is calculated by the CCSD theory. 
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FIG. 2. Potential energy curves (in relative energy) for the ground state of H2, calculated by 
the spin-restricted HF theory and KS-DFT (with various XC functionals). The exact potential 
energy curve is calculated by the CCSD theory. The zeros of energy are set at the respective 
spin-unrestricted dissociation limits. 
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FIG. 4. Same as Fig. [3j but in relative energy. The zeros of energy are set at the respective 
spin-unrestricted dissociation limits. 
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FIG. 5. Occupation numbers of the \a„ orbital for the ground state of H2 as a function of the 



internuclear distance R, calculated by spin-restricted TAO-LDA (with various 9) 
corresponds to spin-restricted KS-LDA. The reference data are the FCI NOONs 
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FIG. 6. Potential energy curves (in relative energy) for the ground state of H2, calculated by the 
spin-restricted (with and without the entropy contributions) and spin-unrestricted TAO-LDA {9 = 
40 mHartree), where the zeros of energy are set at the spin-unrestricted dissociation limit. The 
entropy contributions (in total energy) as a function of the internuclear distance R, calculated by 
spin-restricted TAO-LDA {6 — 40 mHartree), are also shown. 
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FIG. 7. Same as Fig. |6l but for 9 — 7 mHartree. 
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FIG. 8. Potential energy curves (in total energy) for the ground state of N2, calculated by spin- 
restricted TAO-LDA (with various 9). The 8 = case corresponds to spin-restricted KS-LDA. 
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FIG. 9. Same as Fig. [8] but in relative energy. The zeros of energy are set at the respective 
spin-unrestricted dissociation limits. 
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FIG. 10. Occupation numbers of the 3a g orbital for the ground state of N2 as a function of the 
internuclear distance R, calculated by spin-restricted TAO-LDA (with various 9). The 9 = case 
corresponds to spin-restricted KS-LDA. The reference data are the MRCI NOONs (3- 
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FIG. 11. Same as Fig. [TUl but for the 1tt ux orbital. 
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FIG. 12. Potential energy curves (in relative energy) for the ground state of N2, calculated by the 
spin-restricted (with and without the entropy contributions) and spin-unrestricted TAO-LDA {9 — 
40 mHartree), where the zeros of energy are set at the spin-unrestricted dissociation limit. The 
entropy contributions (in total energy) as a function of the internuclear distance R, calculated by 
spin-restricted TAO-LDA {9 = 40 mHartree), are also shown. 
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FIG. 13. Torsion potential energy curves (in relative energy) for the ground state of twisted ethylene 
as a function of the HCCH torsion angle, calculated by spin-restricted TAO-LDA (with various 9). 
The zeros of energy are set at the respective minimum energies. The 9 = case corresponds to 
KS-LDA. 



39 




FIG. 14. Occupation numbers of the tt (lbz) orbital for the ground state of twisted ethylene as 
a function of the HCCH torsion angle, calculated by spin-restricted TAO-LDA (with various 9). 



The 8 = case corresponds to spin-restricted 
NOONs of CASSCF method (HPNO-CAS) 
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FIG. 15. Torsion potential energy curves (in relative energy) for the ground state of twisted ethylene 
as a function of the HCCH torsion angle, calculated by the spin-restricted (with and without the 
entropy contributions) and spin-unrestricted TAO-LDA {9 — 15 mHartree), where the zeros of 
energy are set at the respective minimum energies. The entropy contributions (in total energy) as 
a function of the HCCH torsion angle, calculated by spin-restricted TAO-LDA {0 — 15 mHartree), 
are also shown. 
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FIG. 16. Pentacene, consisting of 5 linearly fuzed benzene rings, is designated as 5-acene. 




FIG. 17. Singlet-triplet energy gap as a function of the acene length, calculated by spin-unrestricted 
KS-DFT (with LDA, BLYP, and B3LYP functionals), using the 6-31G* basis set. The experimental 



data are taken from Refs. 



8fl-l83j, and the DMRG data are taken from Ref. [88 1. 
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FIG. 18. Singlet-triplet energy gap as a function of the acene length, calculated by spin-unrestricted 
TAO-LDA (with various 0), using the 6-31G* basis set. The = case corresponds to spin- 



unrestricted KS-LDA. The experimental data are taken from Refs. 



80|-l83|. 
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FIG. 19. Singlet-triplet energy gap as a function of the acene length, calculated by spin-unrestricted 
TAO-LDA (0 = 7 mHartree), using both the 6-31G* and 6-31G basis sets. 
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FIG. 20. HOMO occupation numbers for the lowest singlet states of ra-acenes as a function of 
the acene length, calculated by spin-restricted TAO-LDA (with various 0)/6-31G*. The 9 = 
case corresponds to spin-restricted KS-LDA. Reference data are the NOONs of the active-space 
variational 2-RDM method 95]. 
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FIG. 21. Active orbital occupation numbers (HOMO-6, HOMO-1, HOMO, LUMO, LUMO+1, 
and LUMO+6) for the lowest singlet states of n-acenes as a function of the acene length, 
calculated by spin-restricted TAO-LDA (9 = 7 mHartree)/6-31G*. 
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TABLES 



TABLE I. Statistical errors (in kcal/mol) of the reaction energies of 30 chemical reactions [20 1 



calculated by TAO-LDA (with various 9 (in mHartree)). The 9 = case corresponds to KS-LDA. 



9 





1 


3 


5 


7 


10 


15 


20 


30 


40 


MSE 


-0.41 


-0.72 


-0.94 


-1.13 


-1.32 


-1.59 


-1.96 


-2.25 


-2.73 


-3.04 


MAE 


8.51 


8.27 


7.75 


7.36 


7.09 


6.95 


7.53 


8.92 


12.28 


15.21 


rms 


11.10 


10.89 


10.31 


9.76 


9.38 


9.16 


9.75 


11.25 


15.20 


19.03 


Max(-) 


-18.31 


-17.43 


-15.65 


-15.73 


-15.92 


-16.55 


-18.63 


-22.61 


-30.65 


-39.72 


Max(+) 


35.68 


35.59 


33.88 


32.18 


30.50 


28.01 


23.95 


20.08 


17.09 


25.45 
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TABLE II. Statistical errors (in A) of EXTS |74j, calculated by TAO-LDA (with various 9 (in 
mHartree)). The 9 = case corresponds to KS-LDA. 
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1 


3 


5 


7 


10 


15 


20 


30 


40 


MSE 


0.004 


0.004 


0.004 


0.005 


0.005 


0.005 


0.006 


0.008 


0.015 


0.030 


MAE 


0.013 


0.013 


0.013 


0.013 


0.013 


0.013 


0.013 


0.014 


0.021 


0.036 


rms 


0.017 


0.017 


0.017 


0.017 


0.017 


0.018 


0.019 


0.021 


0.037 


0.080 


Max(-) 


-0.091 


-0.091 


-0.091 


-0.091 


-0.091 


-0.092 


-0.095 


-0.101 


-0.110 


-0.110 


Max(+) 


0.081 


0.081 


0.080 


0.080 


0.080 


0.080 


0.078 


0.083 


0.222 


0.581 



TABLE III. Singlet-triplet energy gaps (ST gaps) of n-acenes in the polymer limit (re — > oo), 
obtained by nonlinear least-squares fittings of 3 different data sets (20- to 74-acene, 30- to 74- 
acene, and 40- to 74-acene) of the ST gaps calculated by spin-unrestricted TAO-LDA (9 = 7 
mHartree) /6-31G, using a power-law fitting function of the form a + b n~ c . Here, the coefficient of 
determination R 2 is a statistical measure of the goodness-of-fit (R 2 = 1, for a perfect fit). 

Data set ST gap (kcal/mol) a (kcal/mol) b (kcal/mol) c R 2 

20- to 74-acene 0.08 0.077247 ± 0.004245 72.249 ± 0.737 1.1199 ± 0.0038 1.0000 

30- to 74-acene 0.04 0.041992 ± 0.002988 64.485 ± 0.593 1.0806 ± 0.0032 1.0000 

40- to 74-acene 0.03 0.028669 ± 0.001251 61.405 ± 0.284 1.0645 ± 0.0015 1.0000 
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